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Abstract. We give an effective method to compute the entropy for polynomials orthogonal on a segment of the real axis 
that uses as input data only the coefficients of the recurrence relation satisfied by these polynomials. This algorithm is based on 
a series expression for the mutual energy of two probability measures naturally connected with the polynomials. The particular 
case of Gegenbauer polynomials is analyzed in detail. These results are applied also to the computation of the entropy of 
spherical harmonics, important for the study of the entropic uncertainty relations as well as the spatial complexity of physical 
systems in central potentials. 
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1. Introduction. The concept of information entropy in its continuous or discrete form has proved to 
be very fertile in numerous scientific branches because of its flexibility and multiple meanings |17U19ll^Bll38| . 
Indeed, it is used as a measure of disorder in thermodynamics |33| , as a measure of uncertainty in statistical 
mechanics |23 as weu as m classical and quantum information science (2D1 ESI > as a measure of diversity in 
ecological structures, and as a criterion of classification of races and species in population dynamics |23| . 
among others. 

In quantum mechanics, the uncertainty in the localization of a particle in ordinary space is quantitatively 
measured by the so-called position information entropy [Hj 

S p = - J p(r)lnp(r)df, (1.1) 

in a better and more convenient way than the Heisenberg's standard deviation of the quantum-mechanical 
probability density p(r) = \ip(r)\ , where ip(r) is the wavefunction of its dynamical state. Similarly, the 
uncertainty in predicting the momentum of the particle is measured by the momentum information entropy 

2 „ 

of the density ^{p) = ijj(p) , where the Fourier transform ^>{p) of ip(r) is the wavefunction of the same 

state in the dual, conjugate or momentum space. These two quantities describe best |3(JI 131) the extent 
or spread of the position and momentum probability densities, respectively. Moreover, both entropies may 
decrease without bound when the corresponding density becomes more concentrated, i.e. when information 
in the associated space decreases. However, the entropy sum is bounded from below E] 

S p + S~ t > D(l + hi7r), 

where D is the dimensionality of the space (i.e. D = 3 for ordinary space). It expresses the impossibility 
to have a complete information of the position and momentum of the particle simultaneously. This is the 
so-called entropic uncertainty relation, which is a stronger version of the celebrated Heisenberg's uncertainty 
principle, a fundamental law of nature. This fact and the effective implementation of the density functional 
theory of complex many-electron systems j^Z], which uses the single-particle density as the basic variable, 
are responsible for the fact that the study of the entropy has become an ubiquitous tool in some areas (e.g. 
atomic and molecular physics, condensed matter theories). For instance, several maximum entropy methods 
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based on the position and momentum entropies, as well as on their sum, have been developed [161 I21| and 
widely used I3JEDII21 f° r determination of macroscopic quantities of natural systems. Nevertheless, the lack 
of a general theoretical methodology and of accurate numerical algorithms of computation these information 
entropies still prevents this approach from being more widely used. 

The exact determination of the information entropies of complex many-particle systems is a formidable 
task. Only recently a small progress has been achieved using the theory of special functions, which in 
some cases allows to find closed formulas for the information entropies of the simplest 1-dimensional single- 
particle systems and the three-dimensional systems of particles moving in a central or spherically symmetric 
potential. For these systems the wavefunctions are controlled by some classical orthogonal polynomials (such 
as Gegenbauer, Laguerre or Hermite), and the determination of the corresponding information entropies 
boils down naturally to the computation of entropic functionals for sequences of orthogonal polynomials (cf. 
|2s3 EUl ; a state-of-the art of this topic up to 2001 is given in [TT]K 

Namely, for a positive unit Borel measure fx on [—1,1], let 

n 

Pn{x) = 7n J] (* - Cf') , In > , (1.2) 

be the corresponding orthonormal polynomials: 

J Pn{x)p m (x) d(i{x) = S mn , m,neZ + . 
We define the entropy of the polynomials p n (x) as 

E n = E n (/i) = -J pl(x) In (p 2 n (x)) dti(x) . (1.3) 

The presence of these integrals raises two questions. One is the study of their asymptotic behavior when 
n — » oo, which has a special interest in the analysis of the highly-excited (Rydberg) states of numerous 
quantum- mechanical systems of hydrogenic-type |35| . In this sense there have been important contributions 
in the last few years @] El E3 EH EH| > f° r a detailed review, see ^J. A totally different problem is the explicit 
computation of (|1.3fl for every fixed n (up to a certain degree) . Observe that a naive numerical evaluation of 
these functionals by means of quadratures is not convenient: since all the zeros of p n belong to the interval 
of orthogonality, the increasing amount of integrable singularities spoils any attempt to achieve reasonable 
accuracy even for rather small n. 

In this paper we present some theoretical results (Sectional, which allowed us to develop an algorithm 
for an effective and accurate numerical computation of the entropy of polynomials orthogonal on a segment 
of the real axis from the coefficients of the three-term recurrence relation which they satisfy (Section 01 • 
In Section 01 we study in detail the case of Gegenbauer polynomials because of its own interest (as a 
very "representative" class of polynomials l|1.2|l ) and because of their numerous applications; for instance, 
these polynomials control the angular component of the wavefunctions of single-particle systems in central 
potentials (cf. 11 ). In Section|Sl the results of several numerical experiments are discussed, illustrating both 
the accuracy and efficiency of the algorithm proposed here, and comparing it with other computing strategies 
used so far. Finally, the entropy of the spherical harmonics, which measures the spatial complexity of single- 
particle systems and physical systems with central potentials, is computed using the known relationship 
between spherical harmonics and Gegenbauer polynomials (Section 

2. Series representation of the entropy. The entropic functionals l|1.3[) can be restated in terms of 
the logarithmic potential theory. If fi and v are Borel (generally speaking, real signed) measures on C, we 
denote by 



V(z;(i) = \n\z-t\dfi(t) 
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the logarithmic potential of /i, and define the following functionals: 

I\v, /i] — J V(z] v) d/i(z) = — J J In \z — t\ dv(t) d/i(z) , 

is the mutual energy of /u, and u; I(fi) = I[[i,fj] is the logarithmic energy of [i, and 

r i /\ fdv\ , 

R[u, fi\ = — / In — — dv 
J \<W 

is the relative entropy or the Kullback-Leibler information of v and /i. From the Jensen inequality it 
immediately follows that if both /i and v are positive unit measures, then R[v, fi] < 0. 

With the sequence of polynomials 1|1.2|) we can associate naturally two sequences of probability measures 
on [-1,1]: 

1 ™ 

A n = — > 5 An) and dv n (x) — pt(x) d/j,(x) . (2.1) 
n * — ' \f 

j=i 

Both measures are standard objects of study in the analytic theory of orthogonal polynomials. For instance, 
the normalized zero counting measure A„ is closely connected with the n-th root asymptotics of p n , and as 
was shown by Rakhmanov in his pioneering work |28j. v n is associated with the behavior of the ratio p n +\/p n 
as n — > oo. 

With the notations introduced above the entropy (|1.3f) is equivalently rewritten as 

E n = R[u n ,fi\ <0, (2.2) 

or as 

n 

E n = -21n 7 „ + 2^y(C J (n V„) = -21n 7n + 2n7[A„,^ n ]. (2.3) 

A standard procedure for computation of E n employed so far (once the quadratures begin to fail), is 
based on formula l|2.3|l . As it was shown in the potential V(x]v n ) oscillates on [—1,1] around the 
value In 2, which is the Robin (or extremal) constant of this interval. Moreover, the zeros are points 
of local minima for the potential V(x; v n ), hence in order to compute E n (n) we need to sum up the values 
of the logarithmic potential V(z; v n ) at its local minima. Nevertheless, this procedure assumes an explicit 
computation of the zeros of p n , with the drawback of the well-known instability and computational cost of 
this task as n grows large. 

An exception in this sense is an algorithm, proposed in [9] for numerical computation of the entropy 
of Gegenbauer polynomials with integer parameters. It also involves finding zeros of certain polynomials 
generated recursively, but unlike in l|2.3|l . the number of the zeros depends only on the parameter of the 
polynomial and not on its degree. 

In this paper we propose a totally different approach to the computation of the entropic integrals (|1.3fl , 
more in the spirit of standard numerical algorithms for orthogonal polynomials: it uses only the coefficients of 
the recurrence relation satisfied by these polynomials as the input data. It does not involve a solution of any 
nonlinear equation, and it can be carried out by performing matrix multiplication of structured (essentially, 
sparse) matrices. The algorithm is based on the formula contained in Theorem 12. II which seems to be new. 

We denote by L 1 ^ the class of measurable functions on [—1,1], absolutely integrable with respect to \i. 
Let 7fc(x) = cos(fcarccosa;) denote as usual the Chebyshev polynomials of the first kind. With the two 
measures introduced in l|2.1|l we define the double sequence of generalized moments 



Ck,n = / Tk(x)dX n (x), m k , n = / T k (x)du n (x), k,n>0. (2.4) 
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Obviously, \c ktfl \ < 1 and \m k , n | < 1 for all values of k and n. One of the main results of this paper is the 
following 

Theorem 2.1. Assume that fi is a unit Borel measure on [—1,1]. Forn £ N, letp n be the corresponding 
orthonormal polynomial, and the unit measures \ n and v n as defined in h2.1)) . Then for their mutual energy 
the following formula holds: 



I[\ n ,u n ] = ln2 + 2 ]T 



(2.5) 



fc=i 



where the series in the right hand side is convergent. 
Moreover, if we denote 



M n := sup 

xel-i,i] J-i 



pl(x)~pl(t) 



x — t 



dfi(x) < +oo , 



(2.6) 



then, for N <E N we have 



N 



I\X n ,v n ] -ln2-2 



fc=i 



k 



< 



4M„ 
N+ 1 



(2.7) 



Proof. Following we use the Fourier series of the logarithm 18, formula 1.514], 

-ln|l- e -| = £^, 



(2.8) 



k=\ 



valid for almost all tp S [0, n] (see e.g. Theorem 15.2]), which yields a representation of the logarithmic 
kernel 



oo ^ 

\n\x-t\=hx2 + 2^-T k (x)T k (t), 



(2.9) 



fc=i 

where for every 16 [—1,1] the series (in t) converges almost everywhere in [—1,1]. 
On the other hand, by the recurrence relation 

T k+ i(x) = 2xT k (x) - T k -i(x) , T Q (x) = 1, T x {x) = x , 

we have that 

2(t - x)T k (x)T k (t) = q k (x,t) - q k _i(x,t) , q k (x,t) = T k+1 (t)T k (x) - T k+1 (x)T k (t) , 
from where for N £ N, 



(2.10) 



2(t-^)E 



fe=i 



T fc (a;)T fc (t) 
k 



N-l 



-(t- x ) + j2 



Qk(x,t) qN{x,t) 



Hence, 



JY 



2(t - x)J2 TkiX f k(t) 
fc=i 



J k(k + l) N 



<4, t,xe[-l,i] 
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Thus, if / e L 1 we can apply Lebesgue dominated convergence theorem to <|2.9|) in order to assert that 



^/(t)(i- a; )ln|i- : r|^) = /olii2 + 2^^r fc (x), f k = J ^ 



f(t)(t-x)T k (t)d»(t), (2.11) 



and the series in the right hand side is convergent. Furthermore, we can estimate the remainder using that 

T k (x)T k {t) _ q N (x,t) ^ qk(x,t) 

fc=JV+l 



2(t-x) £ 



AT + 1 ' ^ fe(fc + l) 



from where 



so that 



2( t-x) £ ZM^M!) 



< 



2 £ T 7 ^ 



fc=JV+l 

In particular, for 16 [— 1, 1] we may take 



< 



N+l 



|/(*)| £*/*(*) . 



(2.12) 



/(*,■) = 



Pl{t)~P 2 n {x) 
t — x 



and formula (|2.11() yields 



(P 2 n (t) ~ P'W) In |t - a:| 



ln2 T ( | ^( t )_^( a; ))d At (t) + 2f;^ /' (p^(t)-^ (a;) )T fe (t)^(i), 

fe=i 



which using the definitions introduced above can be rewritten as 

Tk(x) 



V(x; u n ) - pl(x)V(x; /x) = (1 - (1)) In 2 + 2 ]T 



fc=i 



Evaluating l|2.13(l at the zeros of p n we obtain 

ncr ) ;^o=in2+25:^f^ 



m k ,n -P n ( x ) I T k (t)dfl(t) 



m k , n , j = 1,2, . . . ,n. 



(2.13) 



k=l 



Summing up these expressions for j = 1, 2, . . . , n, we arrive at l|2.5|) . 
On the other hand, by (|23|l and (|2~T2l . 



2 E 



fc=iV+l 

from where 



(p 2 n (t)-p 2 n (x))T k (t)d»(t) 



< 



N + i y_a 



Pl(t)-Pl(x) 



t — X 



dfj,(t) < 



4M n 
N + l' 



nCf ) ;^)-ln2-2^^i 
fe=i 



2 E 



T k (6 n) ) 



< 



4M» 
7V + 1 



j = 1,2, ... ,71, 
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and the estimate <|2.7[1 follows. □ 

Corollary 2.2. With assumptions of The.ore.m WlS. 



-2 In 



In \ . . V"^ C-k,n m k,n 

2 n ) n ^ k 



(2.14) 



fe=i 



and the error after truncating the series at the N-th term is bounded by the right hand side in [2.7\/ . 

Remark: As our example of Gegenbauer polynomials in Section shows, the bound (|2.7|l is usually too 
pessimistic. 

3. Effective computation of E n . Assume that we have as input data the coefficients of the three-term 
recurrence relation, satisfied by the orthonormal polynomials p n (x) — j n x n + . . . , 



xp n {x) = a n+ ip n+1 (x) + b n p n (x) + a n p n _ 1 (x) . 
with p-i — and po(x) — 1. We form the infinite Jacobi matrix 



(3.1) 



J = 



fb Q oi ...^ 
a± b\ ai 
a 2 b 2 a 3 . . 



V 
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and let J n = J(l : n, 1 : n) denote its principal minor nxn. Here and in what follows we occasionally use the 
MATLAB type notation to refer to elements of a matrix. Furthermore, e„ stands for the infinite (column) 
vector whose n-th element is 1 and the rest is 0, and (a, b) = a H b is the standard scalar product in I 2 or R fc 
(the space where we multiply vectors is always clear from the context). 
The following are very well known facts: 

Proposition 3.1. Letp n be the orthonormal polynomials satisfying the recurrence relation 13.1)) . 

Then, with the notation above, for n > 1, 

i) the zeros cj n \ j = l,...,n, of p n are eigenvalues of J n , and (po(Cj" ) );Pi(Cj" ) ); ■ ■ ■ ;Pri-i(Cj" ) ) T are 
corresponding eigenvectors. In particular, for m £ Nq = N U {0}, 



n 



(n) 



ii) If f is a polynomial then 



Hi) For r > n + 1 and m £ No, 



(e n+ i, /(J)e n+ i) 



= trace of [J n ] 



f(x)pl(x)w(x) dx. 



(r) 



3=1 



m 



(3.2) 



(3.3) 



where A? are the Cotes-C'hristoffel numbers (Gauss quadrature weights) given by 



A 



(r) _ 



-i -i 



(3.4) 
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iv) The leading coefficient of p n satisfies 7„ = (aia 2 ■ ■ ■ a n ) 1 . 

All these facts are classical (see e.g. \32\). Formula <|3.2ll can be found for instance in |^ §4.1.2], where 
it is proved for / £ C[— 1,1]. Identity (13 . 3|) is a straightforward consequence of (|3.4ll . i), and the spectral 
decomposition of the finite sclfadjoint matrix J r . 

For practical computation of the left hand side in 13. 2f) we need the following result: 

Corollary 3.2. If n £ No, m, r £ N and r > n + (m + l)/2, i/ien i/ie (n + l,n+ 1) elements of J m 
and (J r ) m coincide. 

Proof. It is well known that Gauss quadrature formula 



L /(xMx)dx = ^Af/(cj r) ) 

1 3=1 



is exact if the degree of the polynomial / is < 2r — 1. In particular, 

L r 

x m p 2 n (x)w(x)dx = A j r) [Cf] Vl (Cj r) ) 
1 i=i 

is exact if m + 2n < 2r — 1, and the statement follows from H3.2fl - H3.3fl . 1 Q 

Assume that computing the entropy by means of formula l|2.14fi we decide to truncate the series therein 
at k — A £ N. Then as a consequence of the previous corollary, in the right hand side of (|3.2I) we may use 
Tfc(J r ) instead of Tfc(J), with r = n + 1 + [A/2]. Moreover, in the right hand side of If3.2fl it is sufficient 
to know only the (n + l)-th column of f{J r )- Recalling that Chebyshev polynomials satisfy the recurrence 
relation l|2.10fi . we can propose the following algorithm, where /„ stands for the n x n identity matrix: 



Algorithm 1 

(i) Compute In (7„/2 n ) = — In IlJ=i (^ a j) recursively; 
choose N £ N at which truncate the series in ( 12 . 14D ; 
take T (J„) = I n , Ti(J n ) = J n , and iterate 

Tk(Jn) = 2J„ T k ^i(J n ) - r fe _ 2 (J„) , k = 2, . . . , A, 
computing 

(ii) ck, n = tr Tk(J n )/n, k = l,...,N; 

set r = n + 1 + [A/2] ; starting with vq = I r ('-, n + 1) and v\ = J r (:, n + 1) , 
iterate by the recurrence 

v k = 2 J r v k -i - , k = 2, . . . , A, 
computing 

(iii) rrik,n = Vk(n+ 1), k = l, ...,N; 

substitute the results of (i)-(iii) in ( 12 . 141) . terminating the series at k = N. 



Observe that this algorithm starts from the spectral data as the only input, and performs multiplication 
of structured matrices, without solving any kind of equation. This can be efficiently implemented, for 
instance, using the known algorithms for sparse matrix multiplication. 

On the other hand, in order to satisfy conditions of Corollary 13. 21 it is necessary to know in advance the 
truncation term A for the series in (|2.14l) . for which we need an a priori bound on the error. The bound in 
(|2.7I) can be used, but usually it overestimates the error yielding values of A much larger than needed. In 
Section 21 we discuss the selection of the truncation term in the particular case of Gegenbauer polynomials. 



1 This short and elegant prool was suggested by one of the anonymous referees whose contribution we gratefully acknowledge. 
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4. Entropy computation for Gegenbauer polynomials. In this Section we test our approach on 
the computation of the entropy of the Gegenbauer polynomials. 
For A > -1/2 let 

cx= F ( A + 1 ) (41) 

cx v^r(A + i/2)' l4 - ij 

It is easy to verify that w x (x) = c\(l — x 2 )^ 1 / 2 is a positive unit weight on [—1,1]. Let denote the 
Gegenbauer polynomial of degree k and parameter A, orthogonal with respect to w x on this interval, and 
normalized by the value at x = 1, 

Cfr)-(t + *-y, (4-2) 

this is a standard normalization, adopted for instance in pQ and in |32| . Straightforward computation shows 
that 

G x k (x) = ( ^^ + 2^ ) V2 C " (X) = 7feA ^ + l0W6r d6gree temS (4 ' 3) 

are the Gegenbauer polynomials orthonormal with respect to w x (x). 

In this Section we will use the superscript A in the previously introduced notation when we want to 
make the dependence on the parameter A explicit; for instance, 

E* = -j X (Gl{x)f \n(G^x)) 2 w\x)dx. (4.4) 

For the time being, only few explicit formulas for the entropy of orthonormal Gegenbauer polynomials 
are known. Namely, for A = and A = 1 (Chebyshev polynomials of the first and second kind, respectively) 
it is not difficult to prove (cf. ^H3S]) that 

££ = log2-l, E\ = ^— (4.5) 

n + 1 

Furthermore, case A = 2 was studied in [7] and establishing that 

2 f n + 3 \ n 3 - 5?i 2 - 29n - 27 1 / n + 3 \" +2 

n ~° g \3(n + l)J (n + l)(n + 2)(n + 3) n + 2\n+l) ' { ' } 

This formula does not allow to expect a short and elegant expression for all E x , even for integer values of A. 

For A £ N, A > 2, an alternative algorithm has been proposed in P]; it expresses the entropy in terms of 
the zeros of certain polynomials generated recursively, whose number, 2A — 2, depends only on the parameter 
of the polynomial and not on its degree. 

This approach (valid only for integer values of A), is more efficient than the direct computation of E x 
by quadrature when n grows large, and allows also to find constructively the first terms of the asymptotic 
expansion of the entropy when ??. — > oo and A is fixed (cf. |4l I§1 129| ): 

+ £ + tf-l-taj^jj^Lj. (4.7) 

In this Section we will apply the general approach, described in Section to the efficient computation 
of the entropy E x . It is well known that the polynomials G x satisfy the three-term recurrence relation 

x G x (x) = a n +iG^ +1 (x) + a n G^_ 1 (x) , 
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where 



n(n + 2A - 1) 



1/2 



> 0. 



(n + A- l)(n+ A)_ 
In particular, since the leading coefficient 7^ = (a±a 2 ■ ■ . a„) _1 , we have 



-21n( 7 */2") =ln 



j(j + 2A - 1) 
L = l y + A-l)(i + A) 



n 



Consider the values defined in (|2.4|) . By symmetry, 



>k.n 



= for k odd. 



(4.8) 



(4.9) 



Observe that for A > —1/2 the orthogonality weight satisfies the assumptions of Theorem 12. II Thus, taking 
into account l|4.8J) . formula l|2.14|l for the Gegenbauer polynomials reads as 



E n=^ n 



JU + 2A - 1) 
fJi (J+A-I)(j + A) 



2"E 



Col. ^^"1 



- J 2k,n'"'2k,n 



, A>-l/2, neN. 



(4.10) 



fc=i 



PROPOSITION 4.1. For A e N , the series in \4-lU^ is terminating at k = n + A. 

Proof. This is a straightforward consequence of the fact that in this case w x (x)(l — x 2 ) 1 ! 2 is a polynomial 
of degree A, and for k > 2n + 2A we can use orthogonality of Tk in the definition of m& .„ in (|2.4I) . □ 

Taking advantage of (|4.9() . we can use the recurrence formula for even Chebyshev polynomials: To(x) = 1, 
T 2 (x) = 2x 2 - 1, and 

T 2k (x) = 2T 2 (x)T 2k ^ 2 {x)-T 2k - 4 (x), k>2. (4.11) 
Then Algorithm 1 described in Section [3] takes the following form for the Gegenbauer polynomials: 



(i) 



(ii) 



Algorithm 2 
Find -21n(7^/2") by J4~5I) recursively. 

Choose a value A > n + A where to truncate the series in ( 14 . 10D , 

unless A G No ; in such a case, A = n + A. 
Take T (J n ) = I n , J n = T 2 (J n ) = 2J 2 - I„ , and iterate 

T 2 k(J n ) = 2J n T 2 k- 2 (J n ) — T 2 k~4(J n ) , k = 2, . . . , N, 
computing 

c 2fe,„ = tr T 2k (J n )/n, k = l,...,N. 
Set r = n + 1 + [A/2] and J r = T 2 (J r ) = 2 J 2 - I r . 
Starting with vq = I r (:, n + 1) and v 2 = J r (:, n + 1) , 
iterate by the recurrence 

w 2 fc = 2 J r v 2k - 2 - v 2k -i , k = 2, . . . , A, 
computing 



(iii) 



'■2A- 



„ = v 2k (n + 1), k= 1, . 



,A. 



Substitute the results of (i)-(iii) in ( 14 . 10D . 
terminating the series at k = N. 



In the implementation it is convenient to use subroutines for sparse matrix multiplication; as it was 
mentioned, the coefficient in (|4.8|l is better to compute recursively, avoiding possible overflows. 
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In order to obtain an a priori bound for the error (and thus to find the truncation term N) we can 
use the explicit expression for the coefficients rn\ hn . In |35) . were expressed in terms of Wilson 

polynomials of degree n and parameters depending on A and k. Nevertheless, the expression which appears 
there has indeterminacies for integer values of A, which is inconvenient for evaluation. We are interested in 
an alternative formula for m^ k „ • 



When a linearization formula of the type 

T k {x)p n (x) = ^ tj,k,nPj{x) . 



n-\-k 



is available the coefficients mk, n in i|2.4|) can be found observing that Tn^ n — ^n,fc,n- 

Nevertheless, we were 

unable to find the explicit expression in the literature, and we establish a formula for m\ k n based on the 
hypergeometric representation for the Chebyshev and Gegenbauer polynomials. 

Theorem 4.2. For the orthonormal Gegenbauer polynomials G„ the coefficients m^, n defined in \2.4\) 
satisfy 



l 2k,n 



Alternatively, for k > n + X, 



n + A 



3=0 VJ 7 



A(2A + j)„ (-i-A) fe 



J + A (j + A + l) fe ' 



k > 1. 



(4.12) 



^2fc,n = ~ Sin 7rA p > 

j=o 



A (2X + j) n ... , x , . r(A-j-A) 



.7 + A 



^(i + A + 1) 



T(k + j + X + l) 



(4.13) 



In particular, for X € No, m2k,n — /or a/Z fc > n + X + 1. Moreover, for any X > —1/2 and k > n + X, we 
have |m 2 fc+2,n| < \m2k,n\, and 



m 2 fe- 



U2A+1 



(4.14) 



Here (a)jt = T(a + fc)/r(a) denotes the Pochhammer's symbol. 

Proof. The following hypergeometric representation for the Chebyshev and Gegenbauer polynomials is 
well known (see e.g. [22 formulas 8.942 and 8.932]): 



T 2 k(x) = 2*1 
The quadratic transformation 

2 -Pi 

yields 

T 2 k(x) = 2*1 



-2k, 2k 
1/2 



1 - x 



-n, n + 2A 
A + 1/2 



1-x 



-2n, 2n + 2a + 1 
a + 1 



1-z 



— 2-Fl 



-n, n + a + 1/2 
a + 1 



1-z 2 



1/2 



1 - x' 



(2A)„ 



-n/2, n/2 + A 
A + 1/2 



1 - ar 



and Clausen's identity (cf. [21 p. 116, problem 13] or ^] Ch. IV, section 4.3]) 

2 



1*1 



a, b 
a + 6+ 1/2 



3*2 



2a, 26, a + b 
2a + 26, a + 6+ 1/2 
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gives the following representation: 



f(2X) n 



3F2 



[C„(*)] 2 = 

We use these formulas in order to compute the integral 

-C>2fc,n = 



—71, n + 2A, A 
2A, A +1/2 



1 - x' 



^ (-fc)j (k)j ( ^ (-rt)j (n + 2A), (A), f 1 _ , i+i+A _ 1/2 \ 
2- (1/2), i! ^ (2A) i (A + l/2) j i! JJ L X) j 



f (2X)n 

I n! 



2 k n 

EE 



i=0 j=0 

Interchanging the order of summation we get 



(-fc)i (fc)< (-n), (n + 2A), (A), r(i + j + A + l/2) 
(1/2), i\ (2X)j (A + 1/2), j\ r(t + j + A + 1) ' 



/(2A) n \ 2 ^ 



-k, k, j + X + 1/2 
1/2, j + X + 1 



\ (-n)j (n + 2X) j (A),- T(j + A + 1/2) 
(2A), (A + 1/2), j! r(j + A + 1) ■ 



Using the Pfaff-Saalschutz identity, 
the integral Dik,n becomes 

2 n 



a, 6 



c, l + a + 6 — c— A; 



1 = 



(c- a) fc (c- b) k 
(c) fe (c - a - 6) fc 



v 7 j=o 



(1/2 - fc) fc (-j - X) k (-n), (n + 2A), (A), r(j + A + 1/2) 
(l/2) fe {-k - 3 - A) fc (2A), (A + 1/2), j! r(j + A + 1) ' 



Taking into account the normalization factors in (|4.1|l and (|4.3I) we obtain 



x _nl(n + X) T(A + 1) 

m 2fc,n 



T 2k (x)C^(x) 2 (l-x 2 ) x -^ 2 dx 



X (2A)„ V^T(A + l/2) y_! 
(n + A) T(A) (2A)„ (1/2 - fc) fc (-j - A) fe {-n)j (n + 2A), (A), T(j + A + 1/2) 



E 



n! r(A + 1/2) ^ (1/2), (-fc - j - A) fe (2A), (A + 1/2), j\ T(j + A + 1) 



= (-1) 



n+X 



k n + X 



E(-d j 



j=0 



(-j - A) fe , , 1 

(2A + J) n 



3 J (~k-j-X) k 



3 + A 



1! 



E(-D 3 



3=0 



n \ (2X + j) n (-j - X) k 
3 



j + x (j + x + iy 



where we have used standard properties of the Gamma function and Pochhammer's symbol. This proves 
(|4.12|) . Formula l|4.13[) follows easily from H4.12[> and the well known relation r(x)r(l — x) = w/ sin(7ra:) 
applied with x=j + X + l>0. In particular, it shows that |TO2fc+2,n| _• |wi2k,n|> finally, (|4.14() is a 
consequence of (|4. 1 3I> and Stirling formula. □ 

In order to discuss the truncation error in H4.10J1 we need to introduce the following notation: fixed n, 
X > -1/2, A </ N , and N e N, N > n + A, let 



00 A A 
c 2k,n " L 2k,n 



k=N 



12 
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|i?^(iV)| is the absolute error of approximation of if we truncate the series in l|4.10|l after k = N — 1. 
Proposition 14. II shows that it is convenient to take N > n + A. We consider the case A > (for negative A, 
see Remark at the end of this Section): 

Proposition 4.3. Let n 6 N, A > 0, A ^ N, and N e N, N > n + A. Then 



\Rt(N)\ < Fn(N) := 



n(n + A) 
TV 



E 

3=0 



1 



(2A + j)» |(-j-A + l)g_ij 

(n-i)!j! j + A (j + A + 1)tv-i 



(4.15) 



Proof. Let us denote 



n + X (2A + j) n 
j + A (n- j)!jf 



Then by {£12)1 . 



R*(N) = 2nJ2 



l 2fc,n 7— 



k=N 



j=0 fc=A/ 



(— j — A)fe C2fc,n 
(j + A + l)fc ~ 



so that 



\R x n (N)\<2nJ2\A-jY, 

3=0 



k=N 



(j + A + l) fe 



'2fc,n 



k 



<IEi4»iE 



k=N 



(-J - A)a 



(i + A + i) fc 



Taking into account that 

|4„| = (-l)M^ n , and |(-j - A) fc | = (-l)b^- 1 ] (-j - A) fe , for fc > n + [A] + 1, 
where [A] stands for the integer value of A, we get from (|4.17ll . 



' " ( }l " ft J> fro CiTATiW 



But 



A/ 



(-j - A) 



1 (-j — A + 1) 



N-l 



k=0 



(J + A + 1) 



+ 



H - A + 1) 



A/+M 



fc+Af 



2 (j + A + l) W -i 2 (j + A + 1)jv+m 



so that 



\R"n(N)\ < E 4n + = 

TV ^ (j+A+l);v_i nV ; 



(4.16) 



(4.17) 



(4.18) 



and the statement follows. □ 

Given e > we can use (|4.15|) in order to find a (preferably, lowest) value No £ N such that \R^(N)\ < e. 
Obviously, the lower bound for TVo will be TVo = n + [A] + 1 . It is helpful to get also an upper bound for such 
an No, that can be obtained taking advantage of the geometric decay of !F^(N). It is based on the following 



Proposition 4.4. Let n e N, A > 0, A £ N, and N e N, iV > n + A. T/ien /or a// /i e N 0; 



**(JV + /i)< 



(JV + A + h) 2X ' 



(4.19) 
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where T X (N) is defined in \4-15\ ), and 

" V ' ' V 71 V 7 V^ + V ttJV ^ j + X {n-j)\ 3 \ {N-X- 3 ) j {N + X) j 1 7 

is a decreasing function in N , such that 

Proof. We can use the identity 

(-i - A + l) N +h-i i~j ~ X + l)j (A + l)j (-A + l)jv+h-i 



(j + A + l) N +h-i {N + h-X-j)j(N + h + X)j (X + l) N+h -i 



(A^A + ljj (-A + lJ^^-Ajfc 



(JV + ft-A-j)i(JV + fc + A)j (A + l)iv-i (iV + A) h ' 
thus, with notation l|4.1(j[) and by l|4.18f) . 

n(-l)W (—A + 1)jv-i (-AT — X)h sr^r^y ^a (A),(A + 1) 



1 " l jl - JV (A + (iV + A), ^ ij A '> (N + h — X — 3 ) 3 (N + h + A) j 



< n(-l)M (-A+lK-^iV-A)^ ■ A (AM A+ljj 



with 



N (A+1)jv_i (tf + A) h ^ v ■ J,Tl (iV — A — j)j (iV + A)j 

[ N ~ X l h B^N), (4.21) 

-A + 1)jv_i|<A»+A (2A+j)„ (A),(A + 1), 



R A rA n_ » I v A ' -*-7 W— 1 1 \ - 



JV (A + ^ j + A (n - (N-X- j), (N + X), 



r(iV-A)r 2 (A)|Asin(7rA)|^ n + A (2A + j) n (A).,-(A + 1),- 



n 

: iV rfXTV + A) ^ j + A (n - j)! j! (N — X — j)j (N + A), 



where we have used again the identity T(A)r(l — A)sin(7rA) = ir. Alternatively, B X (N) can be represented 
in terms of the following truncating and balanced hypergeometric series, valid for A ^ Z, 

ff A r »n_ 1 l(-A + l)iv-i| n + X (2A) n / -n, n + 2A, A, A 
" l 7 TV (A+1)jv-i (n-1)! A 4 3 ^ 2A,-JV + A + 1,JV + A 

In order to simplify the expression of the error, we may use (cf. j^l P-17]) that for x > y > 1, 

r(») < y^ 1 g 

r(x) ~ a; 31 - 1 e« ' 

Hence, 

(iV — A)ft _ T(JV + A) r(7V — X + h) F(N + X) fN-X + h^'^ 11 ' 1 e 2X 



(N + X)h T(N - X) T(N + X + h) ~ T(N - X) \N + X + hJ (N + X + h) 2X 
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It is straightforward to verify that (x/(x + t)) x_1 is decreasing in x for x,t > 0, so that 

(N-X) h T(N + X) fN-X^-^ 1 e 2X 



< 



\2A 



(N + X) h ~ T(N - X) \N + X ) {N + X + h) 

Gathering this inequality and i|4.22[l in (|4.21|) we obtain the statement of the Proposition. 
□ 

Corollary 4.5. Let e > 0, n € N, A > 0, A £ N, and N € N, N > n + X. Then if for F*(N) defined 

in 



h > max 



-N-X 



(4.23) 



then 



\?Z(N + h)\<e. 



Hence, if we want to find a suitable value of Nq € N for which |i?^(A^)| < s, we can use the following 
procedure: 







Algorithm 3 


(1) 


Take N = n + 


[A] + 1 and compute T*(N); 


(2) 


if T X (N) < e 


, put A^o = N and quit ; 


(3) 


else compute 


h eaual to the r.h.s. of (14 . 23D : 


(4) 


if h > , use 


bisection in order to find the lowest iVo G [N, N + h] PI N such that 




fn(No) < e 





Obviously, we can use a more sophisticated zero-finding method in the procedure above; however, usually 
bisection, which is simple and easy to implement, is sufficient for our needs. 

Remark: A simple alternative for truncation of the series in H4.10(l can be computing E x for different 
integer values of A (say, [A] and [A] + 1) and interpolating the value of for the non-integer A. Nevertheless, 
as numerical experiments in Section [S] show, this approach is not very satisfactory. For small values of A it 
yields large errors, and for large A's the loss in speed computing the entropy at least twice can be compensated 
by larger truncation values N. 

Remark: For small values of the parameter A the error bounds above usually yield large truncation 
values N. This occasionally might justify the use of explicit formulas (|4.12|l and l|4.13|l for computation of 
m 2k n instead of Step 3. As A grows larger, the explicit evaluation of Pochhammer's symbols rapidly becomes 
substantially more time consuming and less accurate than matrix multiplication. 

Furthermore, for — 1/2 < A < the structure of the coefficients rn\ k n yields extremely pessimistic upper 
bounds for the error |i?*(iV)|, with a rate of convergence even lower than 1/N established in Theorem l2.ll 
Nevertheless, in this case the Gegenbauer polynomials are uniformly bounded on [—1,1], and truncation 
error is estimated better using formulas l|2.6|) - (|2.7|) . 

5. Numerical experiments for Gegenbauer polynomials. In this section we discuss briefly the 
performance of the algorithm presented above, and compare it with some alternative algorithms used for 
computing the entropy of Gegenbauer polynomials. 

We will compute E x for several values of A. In all cases Algorithm 2 was implemented in Matlab™ and 
executed on a computer with a single AMD Athlon™ XP2000+ processor, 256 Mb RAM, and running Mat- 
lab 6 under Windows. In particular, specific routines for sparse matrix construction have been used through 
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Degree n 

Fig . 5.1. Error and execution time of the algorithm for computation of . 

Matlab built in functions spdiags and speye. For the experiments no special performance optimization 
techniques have been used, although we implemented vectorization when available. 

One obvious test situation corresponds to A = 0, 1, 2, when the explicit value of the entropy is known (cf. 
formulas 1)4.511 - 1)4.6(1 ). Figure |01 shows that the error of the algorithm (comparing with the exact value 1)4. fijl ) 
is negligible. The execution time, computed as an average of 100 runs of the algorithm, grows geometrically 
with the degree n (due to the proportional growth of the matrices involved) , but still it is below 1 minute 
for n as large as 500. 

It is also illustrating to compare E^ with their asymptotic expansion 1)4.7)1 truncated after the first and 
the second terms (Fig. I5.2|l . 

In order to illustrate the performance of Algorithm 2 we compare it with the following procedures to 
compute £„: 

• By formula 1)4.4)1 using adaptative quadrature implemented in Mathematica™ 4.2; 

• By formula 1)4.4)1 using functions quad and quadl of Matlab 6, applied to explicit expressions of the 
polynomials with coefficients computed using Mathematica 4.2 with exact arithmetics; 

• By the algorithm described in [5]. This approach is valid for integer values of A only. 

In Table l5~D we compare the errors and execution times of the procedures above with those of Algorithm 2. 
The execution time is taken as the average of 100 runs of the corresponding algorithms. 

As it was mentioned above, the bound in 1)4.15(1 usually overestimates the error. For that purpose we 
truncate first the series in ()4.10)) at N such that !F^(N) is not greater than the machine epsilon; the corre- 
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Degree n 

Fig. 5.2. Entropy E^ for A = 2 (marked with '*'), A = 10 (marked with 'o') and A = 20 (marked with 'x') compared with 
its limit £q (top) and the asymptotic expression E$ + E$/n (bottom). 

sponding value of E^ is assumed as the "true" value of the entropy. We compare it with the approximation 
of that we obtain if we truncate the series in i|4.10[l at N given by Algorithm 3 fFig.lQJl. 

Finally, it is tempting to avoid the question of truncation error in H4.10[l by applying Algorithm 2 to 
A e No and computing E^ for non-integer values of the parameter by interpolation. As experiment, we 
compute -E200 f° r half- integer values of A using two strategies. In Figure IQ1 we observe the execution time 
(in seconds) for truncating the series in (|4.1U|) in the way that F^oo — 10~ 6 (dots), and for interpolating 
i?2 00 by cubic splines using the values of the entropy for A ± 1/2, A ± 3/2 (diamonds). Asterisks represent 
the errors of interpolation. As we see, interpolation usually yields large errors for small values of A, where it 
still could be competitive, since for large A's we are penalized by the time invested in computing the entropy 
at least twice. 

6. Computation of the entropy of spherical harmonics. In this section, to show the usefulness 
of our computational algorithm as well as the close connection of the entropy of Gegenbauer polynomials 
analyzed in detail in the three previous sections, we determine the spatial complexity of some quantum- 
mechanical prototype and real systems with central potentials (rigid rotator, harmonic oscillator, hydrogen 
atom, Rydberg atoms, some diatomic molecules, etc.) by means of the entropy of the spherical harmonics, 

Si, m [Y] :=- /|li m (fi)| 2 ln|y ; , m (r!)| 2 dr!, 
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Method 


Absolute error 


Time (sec) 


n= 10 


n = 25 


n = 50 


n = 100 


71 = 10 


n = 25 


n = 50 


n = 100 


(i) 


4.0 x 10~ 7 


1.4 x 10" 4 


5.8 x 10~ 4 


1.4 x 10" 3 


0.30 


0.55 


1.13 


2.51 


(ii) 


4.3 x 10~ 6 


1.3 x 10~ 5 


5.7 x 10~ 5 


1.1 x 10~ 2 


0.08 


0.28 


0.89 


2.61 


(iii) 


7.7 x 10~ 8 


1.4 x 10" 7 


8.6 x 10~ 7 


5.2 x 10~ 6 


0.11 


0.31 


0.74 


3.01 


(iv) 








1.8 x 10~ 26 


1.3 x 10~ 7 


0.006 


0.008 


0.01 


0.02 


(v) 


4.4 x 10~ 15 


2.2 x 10" 16 


6.0 x 10" 15 


2.7 x 10~ 15 


0.002 


0.007 


0.03 


0.25 



Table 5.1 

Absolute error (left half) and time of computation (in seconds) of E% for n = 10, 25, 50, 100 by means of (i) adaptative 
quadrature of Mathematica 4.1 with extended precision (left half) or using floating point arithmetics (right half), (ii) quad 
and (iii) quadl functions of Matlab, (iv) algorithm from implemented in Mathematica 4-1 with extended precision, and (v) 
Algorithm 2. 



where f2 = (0,(p), dfi = smOdddip, with < 9 < it and < <p < 2ir, and Y^ m (£l) denotes the spherical 
harmonics which depend on the orbital and azimuthal quantum numbers, I and m respectively. It is well- 
known that the principal quantum number n € No, together with I € No and m G Z, completely characterize 
a single-particle system with a central potential. Moreover, for a given n the orbital quantum numbers 
I < n — 1, and for a given / the azimuthal quantum number —l<m< +1. 

The spatial or angular wavefunction of the system, which defines its bulky shape, can be expressed in 
terms of Gegenbauer polynomials as (cf. |36j ) 



where are Gegenbauer polynomials normalized as in l|4.2|) . and 



N, 



(; + i)(/-H)![r(|m| + i)]' 

2 i-3|fi»[ 7r 2 (l + \ m \)\ 



1/2 



Then, taking into account relation Q4.3JI . the entropy ^^[F] is expressed in terms of the entropy of the 
Gegenbauer orthonormal polynomials, defined in Q4.4[l . as 



Si, m [Y] =ln 



2tt 



c \m\+l/2 



E] m l + l /2 - \m\ 

I— \m\ I I 



■2r(l + \,n\ +:i)- 1 il + ^-2[u2- T -^— 



where c\ is defined in (|4.1|) . 

Thus, we can apply Algorithm 2 in order to compute the entropy of the spherical harmonics jSz, m [y] for 
different values of the quantum numbers I and m. In Fig. 16.11 values of 6*200, m\X] are computed for integer 
values of < m < 200. This figure illustrates that for a given I the entropy is higher around the center of the 
manifold of azimuthal quantum numbers m = —I, — I + 1, — I + 2, . . . , / — 1, 1, than at its extremes, indicating 
that the spherical harmonics are much more localized for the largest values of |m|. Moreover, the entropy 
is approximately constant in the interval —1/2 < m < 1/2, and then it monotonically decreases when \m\ 
grows up to its largest allowed value I. The origin of this intriguing phenomenon is the delicate interplay of 
the sinus factor and the Gegenbauer polynomial involved in the spherical harmonics, which deserves further 
numerical investigation. 
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N 

Fig. 5.3. Error in computing i?2oo by \4-lU\l truncating the series at N (dots) compared with the error bound F^qq^N) 
(solid line) for A = 1.5 (top) and A = 21.5 (bottom). 
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